#ifndef incl_LagrangeElement_h
#define incl_LagrangeElement_h

#include <vector>
#include "MaterialBase.h"
#include "InternalVariables.h"
#include "ElementTypeData.h"
#include "SolutionDataSolid.h"
#include "GeomDataLagrange.h"
#include "headersEigen.h"
#include "petscsystypes.h"

using std::vector;
using std::cout;
using std::endl;


using Eigen::VectorXd;
using Eigen::MatrixXd;


class ElementBase
{
  private:
    static int  elemcount;

  public:

    //member variables

    ElementShape  ELEM_SHAPE;
    int  ELEM_TYPE;
    string  ELEM_FORMULATION;

    int  elmType, matType, secType, subdomId, elmTypeNameNum;
    int  matId, npElem, ndim;
    int  ndof, nsize, nivGP, nGP, elenum, nlbfP;

    double *elmDat, *matDat;
    double  elemError, elemVol, elemVolCur, elemVolOrig;
    double  pres, charlen, dtCrit;

    vector<double>  resi, primvar, vals2project, resi2;
    InternalVariables  ivar;

    vector<int>  nodeNums, nodeNumsPres, forAssyVecPres, globalDOFnums, nodeInLocalElem;
    vector<PetscInt> forAssyVec;

    SolutionDataSolid  *SolnData;
    GeomDataLagrange   *GeomData;
    MaterialBase       *MatlData;
    ElementTypeData    *ElemTypeData;

    // entries for mixed and Bbar formulations
    VectorXd  Nc, dNc_dx, dNc_dy, dNc_dz, Rp;
    VectorXd  presDOF, presDOFprev, presDOFprev2;
    MatrixXd  Kup, Kpu, Kpp, Kut, Ktu;

    //member functions

    ElementBase();

    virtual ~ElementBase();

    int getDimension()
    { return ndim; }

    void  setSubdomainId(int sid)
    {  subdomId = sid; return;  }

    int getSubdomainId()
    {  return  subdomId;  }

    ElementShape getElementShape()
    {  return  ELEM_SHAPE;  }

    int getNodesPerElement()
    {  return  npElem; }

    int getNdofPerNode()
    { return ndof; }

    int  getNdofPerElement()
    {  return  nsize;  }

    int  getElementTypeNumber()
    {  return  elmType;  }

    int  getMaterialTypeNumber()
    {  return  matType;  }

    int  getSectionTypeNumber()
    {  return  secType;  }

    std::vector<int>&  getNodeNumbers()
    {  return  nodeNums; }

    std::vector<PetscInt>&  getVectorForAssembly()
    {  return  forAssyVec; }

    virtual ElemTypeNum getElmTypeNameNum()
    {  return ELEM_DUMMY; }

    virtual void printStiffnessMatrix();

    virtual void printForceVector();

    virtual void prepareElemData();

    virtual void prepareElemData2()
    { cout << "   'prepareElemData2()' is not defined for this element!\n\n"; return; }

    virtual void prepareElemData(vector<myPoint>& nodeCoords)
    { cout << "   'prepareElemData' is not defined for this element!\n\n"; return; }

    virtual void printPrimVariable();

    double getError()
    { return  elemError;   }

    double getVolume()
    { return  elemVol;  }

    virtual void initialiseDOFvalues()
    { cout << "   'initialiseDOFvalues' is not defined for this element!\n\n"; return; }

    virtual int calcOutput(double u1, double v1)
    { cout << "   'calcOutput' is not defined for this element!\n\n"; return 0; }

    virtual void setnivGP();

    virtual void initialiseIntVar();

    virtual int calcInternalForces()
    { cout << "   'calcAndAssyIntForceVec' is not defined for this element!\n\n"; return 0; }

    virtual int calcLoadVector(VectorXd& Flocal)
    { cout << "   'calcLoadVector' is not defined for this element!\n\n"; return 0; }

    virtual int  calcStiffnessAndResidual(MatrixXd& Klocal, VectorXd& Flocal)
    { cout << "   'calcStiffnessAndResidual' is not defined for this element!\n\n"; return 0; }

    virtual int  calcStiffnessAndResidualMixed(MatrixXd& Kuu, MatrixXd& Kup, MatrixXd& Kpu, MatrixXd& Kpp, VectorXd& Flocal1, VectorXd& Flocal2)
    { cout << "   'calcStiffnessAndResidualMixed' is not defined for this element!\n\n"; return 0; }

    int calcStiffnessAndResidual_Formulation_Displacement_2D(MatrixXd& Klocal, VectorXd& Flocal);

    int calcStiffnessAndResidual_Formulation_Displacement_3D(MatrixXd& Klocal, VectorXd& Flocal);

    virtual int  calcMassMatrix(MatrixXd& Mlocal, bool MassLumping)
    { cout << "   'calcMassMatrix' is not defined for this element!\n\n"; return 0; }

    virtual int  calcMassMatrixPressure(MatrixXd& Mlocal, bool MassLumping)
    { cout << "   'calcMassMatrixPressure' is not defined for this element!\n\n"; return 0; }

    virtual int  calcResidual(VectorXd& Flocal)
    { cout << "   'calcResidual' is not defined for this element!\n\n"; return 0; }

    virtual int  calcResidualPressure(VectorXd& Flocal)
    { cout << "   'calcResidualPressure' is not defined for this element!\n\n"; return 0; }

    virtual double  calcCriticalTimeStep(bool flag)
    { cout << "   'calcCriticalTimeStep' is not defined for this element!\n\n"; return 0; }

    virtual int  solveForPressure()
    { cout << "   'solveForPressure' is not defined for this element!\n\n"; return 0; }

    virtual int  applyDirichletBCs(MatrixXd& Klocal, VectorXd& Flocal);

    virtual int  applyDirichletBCsMixed(int offset, MatrixXd& Kuu, MatrixXd& Kup, MatrixXd& Kpu, MatrixXd& Kpp, VectorXd& Flocal1, VectorXd& Flocal2);

    virtual int  applyDirichletBCs2field(int var3_offset, MatrixXd& K11, MatrixXd& K12, MatrixXd& K21, MatrixXd& K22, VectorXd& Flocal1, VectorXd& Flocal2, vector<PetscInt>& forAssyVar1, vector<int>& forAssyVar2, VectorXd& var1applied, VectorXd& var2applied);

    virtual void contourplot(int, int, double, double)
    { cout << "   'contourPlot' is not defined for this element!\n\n"; return; }

    virtual void elementContourplot(int, int, int)
    { cout << "   'elementContourplot' is not defined for this element!\n\n"; return; }

    virtual double computeVolume(bool init = false)
    { cout << "   'computeVolume' is not defined for this element!\n\n"; return 0.; }

    virtual void plotGaussPoints(int, bool defFlg = false)
    { cout << "  'plotGaussPoints' is not available for this element!\n\n"; return; }

    virtual void projectToNodes(bool, int, int, int)
    { cout << "  'projectToNodes' is not available for this element!\n\n"; return; }

    virtual void projectInternalVariable(bool, int, int, int, double*)
    { cout << "  'projectInternalVariable' is not available for this element!\n\n"; return; }

    virtual void projectStrain(bool, int, int, int, double*)
    { cout << "  'projectStrain' is not available for this element!\n\n"; return; }

    virtual void projectStress(bool, int, int, int, double*)
    { cout << "  'projectStress' is not available for this element!\n\n"; return; }

    virtual  void computeEnergy(int, int, VectorXd& )
    { cout << "  'computeEnergy' is not available for this element!\n\n"; return; }

    virtual  int calcError(int index)
    { cout << "  'calcError' is not available for this element!\n\n"; return -1; }

    void  computeDefGrad2D(VectorXd& dN_dx, VectorXd& dN_dy, MatrixXd& F);
    void  computeDefGrad2DPrev(VectorXd& dN_dx, VectorXd& dN_dy, MatrixXd& F);
    void  computeDefGrad2DCur(VectorXd& dN_dx, VectorXd& dN_dy, MatrixXd& F);

    void  computeDefGrad(VectorXd& dN_dx, VectorXd& dN_dy, VectorXd& dN_dz, MatrixXd& F);
    void  computeDefGradPrev(VectorXd& dN_dx, VectorXd& dN_dy, VectorXd& dN_dz, MatrixXd& F);
    void  computeDefGradCur(VectorXd& dN_dx, VectorXd& dN_dy, VectorXd& dN_dz, MatrixXd& F);
};

#endif //incl_Lagrange_Element_h

